The role of fixed delays in neuronal rate models 
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Fixed delays in neuronal interactions arise through synaptic and dendritic processing. Previous 
work has shown that such delays, which play an important role in shaping the dynamics of networks 
of large numbers of spiking neurons with continuous synaptic kinetics, can be taken into account 
with a rate model through the addition of an explicit, fixed delay. Here we extend this work to 
account for arbitrary symmetric patterns of synaptic connectivity and generic nonlinear transfer 
functions. Specifically, we conduct a weakly nonlinear analysis of the dynamical states arising via 
primary instabilities of the stationary uniform state. In this way we determine analytically how 
the nature and stability of these states depend on the choice of transfer function and connectivity. 
While this dependence is, in general, nontrivial, we make use of the smallness of the ratio in the 
delay in neuronal interactions to the effective time constant of integration to arrive at two general 
observations of physiological relevance. These are: 1 - fast oscillations are always supercritical for 
realistic transfer functions. 2 - Traveling waves are preferred over standing waves given plausible 
patterns of local connectivity. 

PACS numbers: 87.19.lj, 87.19.lp, 05.45.-a, 84.35.+i, 89.75.-k 



I. INTRODUCTION 

When studying the collective dynamics of cortical neu- 
rons computationally, networks of large numbers of spik- 
ing neurons have naturally been the benchmark model. 
Network models incorporate the most fundamental phys- 
iological properties of neurons: sub-threshold voltage dy- 
namics, spiking (via spike generation dynamics or a fixed 
threshold), and discontinuous synaptic interactions. For 
this reason, networks of spiking neurons are considered to 
be biologically realistic. However, with few exceptions, 
e.g. [31 m |S], network models of spiking neurons are not 
amenable to analytical work and thus constitute above 
all a computational tool. Rather, researchers use reduced 
or simplified models which describe some measure of the 
mean activity in a population of cells, oftentimes taken 
as the firing rate. Firing-rate models are simple, phe- 
nomenological models of neuronal activity, generally in 
the form of continuous, first-order ordinary differential 
equations [5]. Such firing-rate models can be analyzed 
using standard techniques for differential equations, al- 
lowing one to understand the qualitative dependence of 
the dynamics on parameters. Nonetheless, firing-rate 
models do not represent, in general, proper mathematical 
reductions of the original network dynamics but rather 
are heuristic, but see [7\ As such, there is in general 
no clear relationship between the parameters in the rate 
model and those in the full network of spiking neurons, 
although for at least some specific cases quasi- analytical 
approaches may be of value [8]. It therefore behooves 
the researcher to study rate models in conjunction with 



network simulations in order to ensure there is good qual- 
itative agreement between the two. 

Luckily, rate models have proven remarkably accu- 
rate in capturing the main types of qualitative dynam- 
ical states seen in networks of large numbers of asyn- 
chronously spiking neurons. For example, it was shown 
in [1, 2J that the addition of an explicit delay in a rate 
equation was sufficient to describe the emergence of fast 
oscillations prevalent in networks of spiking neurons with 
dominant inhibition. In such networks, effective delays 
due to the continuous synaptic kinetics can lead to os- 
cillations even in the absence of any explicit delay in 
the neuronal interactions. When the pattern of synaptic 
connectivity depends on the distance between neurons, 
then this effective delay can also lead to the emergence 
of waves. This is certainly a relevant case for local cir- 
cuits in cortical tissue, where the likelihood of finding a 
connection between any two neurons decreases as a func- 
tion of the distance between them, e.g. [S]. In [TJ [2], 
the authors studied a rate model with fixed delay on a 
ring geometry with two simplifying assumptions. First 
they assumed that the strength of connection between 
neurons could be expressed as a constant plus the co- 
sine of the distance between the neurons. Secondly, they 
assumed a linear rectified form for the transfer function 
which relates inputs to outputs. These assumptions al- 
lowed them to construct a detailed phase diagram of dy- 
namical states, to a large degree analytically. In addition 
to the stationary bump state (SB) which had been stud- 
ied previously [lOllIT], the presence of a delay led to two 
new states arising from primary instabilities of the sta- 
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tionary uniform state (SU): an oscillatory uniform state 
(OU) and a traveling wave state (TW). Secondary bi- 
furcations of these three states (SB,OU,TW) led to yet 
more complex states including standing waves (SW) and 
oscillatory bump states (OB). Several regions of bistabil- 
ity between primary and secondary states were found, in- 
cluding OU-TW, OU-SB and OU-OB. They subsequently 
confirmed these results through simulations of networks 
of Hodgkin-Huxley neurons. Despite the good agreement 
between the rate equation and network simulations, two 
main issues remain unresolved. 

• The rate equation predicted that the primary insta- 
bility of the SU state to waves should be to traveling 
waves, while in the network simulations standing 
waves were robustly observed. 

• The linear-threshold transfer function, albeit 
amenable to analysis, nonetheless leads to degen- 
erate behavior at a bifurcation point. Specifically, 
any perturbations with a positive linear growth rate 
will continue to grow until the lower threshold of 
the transfer function is reached. This means that 
the amplitude of new solution branches at a bi- 
furcation is always finite, although the solution it- 
self may not be subcritical. In a practical sense 
then, this means that it is not possible to assess 
whether a particular solution, for example oscilla- 
tions or bumps, will emerge continuously from the 
SU state as a parameter is changed, or if it will ap- 
pear at finite amplitude and therefore be bistable 
with the SU state over some range. 

In order to address these issues, and provide a more com- 
plete analysis of the role of fixed delays in neuronal tissue, 
we here study a rate equation with delay without impos- 
ing any restrictions on the form of the transfer function 
beyond smoothness or on the shape of the connectiv- 
ity kernel beyond being symmetric. What we show is 
that the nature and stability of the solutions arising via 
primary instabilities of the SU state (oscillations, bumps 
and waves) depend on nonlinear combinations of the first 
three derivatives of the transfer function and the first 
three spatial Fourier coefficients of the connectivity ker- 
nel. Thus while the presence of a fixed delay alone is 
sufficient to generate oscillations and waves, whether the 
oscillations are bistable with the unpatterned state and 
which type of waves (standing waves or traveling waves) , 
appear as a stable solution, depend crucially on the trans- 
fer function and pattern of synaptic connectivity. We will 
discuss this in great depth in the results section. 

We would like to emphasize the fact that we are in- 
terested in the effect of a fixed delay on the dynam- 
ics of a local patch of cortical tissue. In fact, de- 
lays in the nervous system are most often associated 
with transmission delays, i.e. delays due to the fi- 
nite velocity propagation of action potentials along ax- 
ons. Indeed, this type of propagation delay, which de- 
pends linearly on the distance between any two neu- 
rons, has been the topic of much theoretical study, e.g. 



[ni[Tl[Tl[ISl[TSJ[I71[THJ[Til2n]. Localized solutions of 
integro-differential equations describing neuronal activ- 
ity, including fronts and pulses, are affected by distance- 
dependent axonal delays [HJ [T5J [T71 Specifically, 
the velocity of propagation of the localized solution is 
proportional to the conduction velocity along the axon 
for small conduction velocities, while for large conduc- 
tion velocities it is essentially constant. This reflects 
the fact that the propagation of activity in neuronal tis- 
sue is driven by local integration in which synaptic and 
membrane time constants provide the bottleneck. Also, 
allowing for different conduction velocities for separate 
excitatory and inhibitory populations can lead to bifur- 
cations of localized bump states to breathers and trav- 
eling pulses [I8|. Pattern-forming instabilities have also 
been studied in the presence of distance-dependent de- 
lays. The presence of propagation delays in rate models 
can lead to oscillations and waves [12l [13] . The weakly 
nonlinear dynamics of waves in spatially extended rate 
models, i.e. describing large-scale (on the order of cen- 
timeters) activity, is described by the coupled mean-field 
Ginzburg-Landau equations |14j . and thus exhibits the 
phenomenology of small amplitude waves familiar from 
other pattern forming systems [21] . 

On the other hand, fixed delays, and their effect on the 
dynamics of large numbers of recurrently coupled neu- 
rons, have received relatively little theoretical attention. 
Some exceptions include work on the role of global feed- 
back delay in shaping the power spectrum of a network 
driven by noisy inputs [22] , and the effect of distributed 
delays in a mean-field corticothalamic model [23|. More 
akin to the work we present here is the recent study of 
the effect of two distinct delays in a spatially homoge- 
neous Wilson-Cowan model [23]. There the authors were 
able to compute oscillatory solutions analytically given a 
Heaviside transfer function of the rate variables. 

More generally, fixed delays, which are likely due to 
synaptic and dendritic integration, and conduction de- 
lays due to the propagation of action potentials along 
the axon, are both present in real neuronal systems. Im- 
portantly, this means that the delay in neuronal interac- 
tions at zero distance is not zero. In fact, fixed delays 
are always observed in paired intracellular recordings in 
cortical slices. The latency from the start of the fast 
rising phase of the action potential to the start of the 
post-synaptic current (or potential) has been measured 
for pairs of pyramidal cells in rat layers 3 to 6 and is 
on the order of milliseconds, see [25] for a recent review. 
Recordings from cat cortex and between pyramidal cells 
and other cells including spiny cells and interneurons in 
the rat cortex also reveal fixed delays which are rarely 
less than a millisecond. These delays are seen when neu- 
rons are spatially adjacent, indicating that axonal prop- 
agation is not an important contributing factor. On the 
other hand the speed of propagation of action potentials 
along unmyelinated axons in mammals is on the order of 
10""'^ — 10^ ni/s, which means a delay of 0.1-10 ms for neu- 
rons separated by 1 millimeter [2S1 127| . Thus fixed delays 



3 



and conduction delays are of similar magnitude within 
a local patch of cortex and both would be expected to 
shape the dynamics of non-steady activity, i.e. neither is 
negligible. Here we have decided to focus on fixed delays, 
as in previous work [TJ [2] , due both to their physiological 
relevance and prevalence in networks of spiking neurons. 

Thus in what follows we will study a rate equation 
with fixed delay and spatially modulated connectivity. 
In section |ll] we formulate the model and conduct a lin- 
ear stability analysis of the SU state. In section [TTT| we 
conduct a weakly nonlinear analysis for the four possi- 
ble primary instabilities of the SU state (asynchronous 
unpatterned state in a network model), thereby deriv- 
ing amplitude equations for a steady, Turing (bumps), 
Hopf (global oscillations), and Turing-Hopf (waves) bi- 
furcations. We will focus on the delay-driven instabili- 



ties, i.e. Hopf and Turing-Hopf. Finally, in section IV we 



will study the interactions of pairs of solutions: bumps 
and global oscillations and global oscillations and waves 
respectively. 



II. THE MODEL 

We study a heuristic equation describing the activity 
of a small patch of neural tissue consisting of two pop- 
ulations of recurrently coupled excitatory and inhibitory 
neurons respectively. Our formulation is equivalent to 
the Wilson-Cowan equations without refractory period 
[6], and with spatially dependent synaptic connectivity 
which was studied originally in PS^. Additionally, we 
consider a fixed delay in the neuronal interactions. Given 
these assumptions, the full Wilson-Cowan equations are 



eT-e = -Te + ^ei / dyJee{\x-y\)re{y,t-de)~ / dy Jei{\x - y\)ri{y , t - di) + 1^ ] , 



(1) 



Tih = + / dyJie{\x - y\)re{y,t - de) - / dyJu{\x-y\)n{y,t~di) + I, 



(2) 



In the original formulation j6j, re{x,t) and ri{x,t) repre- 
sent the average number of active cells in the excitatory 
and inhibitory populations respectively, in this case at a 
position x and at a time t. The time constant Te (r^) is 
roughly the time it takes for a an excitatory (inhibitory) 
cell receiving "at least threshold excitation" [S] to gener- 
ate a spike. This can reasonably be taken as the mem- 
brane time constant which is generally on the order of 10- 
20 ms. The functions ^a{x){a = e, i) are usually taken to 
be sigmoidal. Specifically, if all neurons in the population 
receive equal excitatory drive, and there is heterogeneity 
in some parameter across neurons, e.g. the threshold to 
spiking, which obeys a unimodal distribution, then the 
fraction of active neurons is just the integral over the 
distribution, up to the given level of excitation. The in- 
tegral of a unimodal distribution is sigmoidal. In Eqs[T]- 
[2] the functions Jab{\x\){a = e,i){b = e,i) represent the 
strength of synaptic connection from a neuron in popula- 
tion 6 to a neuron in population a separated by a distance 
X. Here the neurons are arranged in one dimension on 
a domain f2. Input from excitatory (inhibitory) cells is 
furthermore delayed by a fixed amount de (di), which, 
as we have discussed in the introduction, is on the order 
of one millisecond. Finally, the excitatory and inhibitory 
populations are subject to an external drive of strength 
le and li respectively. 

A general analysis of Eqs{l][2] would be technically ar- 
duous although it is a natural extension of the work pre- 
sented here. Rather, we choose to study the dynamics 



of this system under the simplifying assumption that the 
excitatory and inhibitory neurons follow the same dy- 
namics, i.e. Tg = Ti = T, de = di = d, J^e = Jie = Je, 

Jei = Jiz = Jt, *e = ^ <i>, le ^ It ^ I- If this the 
case, then rg = r; = r and the variable r follows the 
dynamics given by 



r{x,t) = -r(a;,t)-|-<i>l 



27r 



dyJi\x-y\)riy,t-D)+I 



(3) 

where we have chosen the domain to be a ring of nor- 
malized length L = 2tt. Furthermore, we have re-scaled 
time by the time constant r. The normalized delay is 
therefore D = d/r, which is the ratio of the effective delay 
in neuronal interactions to the integration time constant 
and should be much less than one in general. The synap- 
tic connectivity expressed in terms of the excitatory and 
inhibitory contributions is J(|a;|) = Je(|a;|) — ^i(|a:|) and 
thus represents an effective mixed coupling which may 
have both positive and negative regions. Eq(3] with the 
choice of <&(/) = / for x > and otherwise and with 
J{x) = Jq + Ji cos (x) is precisely the model studied in 
[H 12] . We now wish to study Eqjs] for arbitrary choices 
of «>(J) and J{x). 

In presenting EqjSjwe have relied on the heuristic phys- 
iological motivation first put forth in jGj. Nonetheless, as 
a phenomenological model, the terms and parameters in 
Eq|3]may have alternative and equally plausible interpre- 
tations. Indeed, the variable r is often thought of as the 
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firing rate as opposed to the fraction of active cells, in 
which case the function can be thought of as the 

transfer function or fl curve of a cell. Experimentally 
this function has been found to be well approximated by 
a power-law nonlinearity with a power greater than one 
[29l [30] . Modeling studies show that the same nonlinear- 
ity applies to integrate-and-fire neurons and conductance 
based neurons driven by noisy inputs |30j . Therefore it 
may be that such a choice of $ leads to better agreement 
of Eq|3]with networks of spiking neurons and hence with 
actual neuronal activity. More fundamentally, we may 
ask if choosing $ as a sigmoid or a power law qualita- 
tively alters the dynamical states arising in Eq|3] This is 
precisely why we choose here not to impose restrictions 
on $ but rather conduct an analysis valid for any <i>. How 
the choice of $ affects the generation of oscillations and 
waves is an issue we will return to in the corresponding 
sections of this paper. 



A. Linear stability analysis 

Stationary uniform solutions (SU) of Eqjsjare given by 

i?=$(joi? + /), (4) 

where i? is a constant non-zero rate, Jq is the zeroth order 
spatial Fourier coefficient of the symmetric connectivity 
which can be expressed as 



J{x) = Jo + ^ Jke''''' + c.c. 



(5) 



, k=l 



and k is an integer. Depending on the form of Eq|4] 
may admit one or several solutions. 

We study the linear stability of the SU state with the 
ansatz 



'•(a;,f) = i?+^frfee*'="+" 



(k)t 



(6) 



k=0 



where Srk ^ 1 and the spatial wavenumber k is an in- 
teger due to the periodic boundary conditions. Plugging 
Eq|6]into Eq|3]leads to an equation for the complex eigen- 
value a{k) 



a{k) = -1 + $ Jfee 



-a{k)D 



(7) 



where the slope <i> is evaluated at the fixed point given 
by Eq|4j The real and imaginary parts of the eigenvalue 
a{k) — A(fc) -|- iuj(k) represent the linear growth rate and 
frequency of perturbations with spatial wavenumber k 
respectively. At the bifurcation of a single mode, the 
growth rate will reach zero at exactly one point and be 
negative elsewhere. That is, \{kcr) — for the critical 
mode kcr- Given this, EqjT] yields the dispersion relation 
for the frequency of oscillation of the critical mode 



From Eq|8] it is clear that the wavelength of the criti- 
cal mode depends crucially on the synaptic connectivity. 
In particular, the spatial Fourier coefficients of the con- 
nectivity kernel J{x) depend on the wavenumber k, i.e. 
Jk — J{k). Thus, the critical wavenumber is, in effect, 
selected by the choice of connectivity kernel. It is in this 
way that the nature of the instability depends on the 
synaptic connectivity at the linear level. 

Depending on the values of uj and kcr in EqjS] at the 
bifurcation from the SU state, four types of instabilities 
are possible: 

• Steady (w = 0, kcr = 0): the instability leads to a 
global increase in activity. 

• Turing {lj = 0, kcr 7^ 0): the instability leads to a 
stationary bump state. 

• Hopf {lo ^ 0, kcr = 0): the instability leads to an 
oscillatory uniform state. 

• Turing-Hopf {uj ^ 0,kcr ^ 0): the instability leads 
to waves. 



For the non-oscillatory instabilities (i.e. a; = 0), Eqj8| 
gives the critical value 



Jk - 1/*' 



(9) 



while for the oscillatory ones Eq[8] is equivalent to the 
system of two transcendental equations 

cZ> = — tantZiZ), (10a) 



-<i> Jk sinujD. 



(10b) 

, = Jk 



Note that we have defined the critical values as Jfc^ 
and ujcr = 

As fixed delays are on the order of a few milliseconds 
and the integration time constant is about an order of 
magnitude larger, the limit of I? —> is a relevant one 
physiologically. This allows us to gain some intuition 
regarding the effect of fixed delays on the dynamics by 
deriving asymptotic results in the limit of small delay. 
Therefore throughout this work we will present asymp- 
totic results, and compare them to the full analytical 
formulas, as well as numerical simulations. In the limit 
of small delay D — > 0, it can be easily shown that an 
instability of the fc*'' spatial Fourier mode occurs at the 
critical value of the coupling 

TT 

with a frequency 



Jk 



(11) 



2D 



(12) 



iuj(kcr) 



-1 $ JkC 



-iuj{kcT)D 



(8) 



FigjTjshows the critical frequency and coupling as a func- 
tion of the delay, up to a delay of 1. The solution obtained 
from the dispersion relation Eqs |10a| and |10b| are given 
by solid lines, while the expressions obtained in the small 
delay limit are given by dotted lines. Thus the expres- 
sions in the small delay limit agree quite well with the 
full expressions even for D = 1. 
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FIG. 1: Top: The critical frequency of oscillatory instabilities 
as a function of the delay D from the dispersion equation 
Eq 10a (solid line) and in the small delay limit (dotted line). 
Bottom: The critical coupling as a function of the delay D 
from Eq 10b (solid) and in the small delay limit (dotted). 



B. An illustrative Phase Diagram 

Throughout the analysis which follows we will illus- 
trate our results with a phase diagram of dynamical 
states. Specifically, we will follow the analysis in [1] |2] 
in constructing a phase diagram of dynamical states as a 
function of Jq and Ji, the first two Fourier coefficients 
of the synaptic connectivity. We will set the higher 
order coefficients to zero for this particular phase dia- 
gram, although we will discuss the effect of additional 
modes in the text. Furthermore, unless otherwise noted, 
for simulations we choose a sigmoidal transfer function 
^(I) = with a = 1.5 and /? = 3. As we vary 

the connectivity in the phase diagram, we also vary the 
constant input / in order to maintain the same level of 
mean activity, i.e. we keep R — 0.1 fixed. For the val- 
ues of the parameters we have chosen here this results in 

/ O.lJo - 0.88. We also take D = 0.1 unless noted 

otherwise 

The primary instability lines for the SU state can be 
seen in the phase diagram, Fig[2] The region in ( Jq, Ji) 
space where the SU state is stable is shown in gray, while 
the primary instabilities, listed above, are shown as red 
lines (color online). In Section III we will provide a 



detailed analysis of the bumps, global oscillations and 
waves (SB, OU and SW/TW) which arise due to the 
Turing, Hopf and Turing-Hopf instabilities respectively. 
The derivation of the amplitude equations is given in Ap- 
pendix |A] as well as a brief discussion of the steady, tran- 
scritical bifurcation which occurs for strong excitatory 
coupling and is not of primary interest for this study. 
Finally, in Section |IV| we will analyze the codimension 
2 bifurcations: Hopf and Turing-Hopf (OU and waves), 
and Turing and Hopf (SU and OU). This analysis will al- 
low us to understand the dynamical states which appear 
near the upper and lower left hand corners of the grey 
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FIG. 2: (Color online) Phase diagram of the rate model Eqjs] 
In each region, the type of solution seen in numerical sim- 
ulations is indicated by a letter code: SU - stationary uni- 
form (grey region), HA - high activity, SB - stationary bump, 
OB - oscillatory bump, SW - standing waves, TW - traveling 
waves. Solid lines indicate analytical expressions. In particu- 
lar, the four possible instabilities of the SU state are depicted 
in red (thick lines correspond to subcritical bifurcations) and 
are given by the linear stability criteria Eqs |9|10"bl The four 
lines emanating from the upper and lower left corners of the 
SU region were determined from a weakly nonlinear analysis 
at the two corners (codimension two points), see section IV. 
The region marked OB corresponds to a mixed mode solu- 
tion of SB-OU, while in the lower left-hand region the OU 
and SW solutions are bistable. Parameters: <l>(x) = , , " „^ 

^ ' l + e 

where a = 1.5 and /3 = 3. We consider the coupling function 
J{x) — Jo + 2Jicosa::. The time delay is -D = 0.1 and the 
input current / is varied so as to keep the uniform stationary 
solution fixed at R — 0.1. 



shaded region in Fig(2) i.e. the SW/OU and OB states. 

III. BIFURCATIONS OF CODIMENSION 1 

As we are interested in creating a phase diagram as 
a function of the connectivity, we will take changes in 
the connectivity as the bifurcation parameter. The small 
parameter e is therefore defined by the expansion 



Jk = Jk + e^A Jfe 



(13) 



The perturbative method we apply, which makes use of 
this small parameter, is called the multiple-scales method 
and is a standard approach for determining the weakly 
nonlinear behavior of pattern-forming instabilities, see 
[2T| . We choose the particular scaling of in the fore- 
knowledge that if the amplitudes of the patterns of in- 
terest are scaled as e, a solvability condition will arise 
at order e^. This solvability condition yields a dynamical 
equation governing the temporal evolution of the pattern. 
See Appendix A for details. Without loss of generality 
we will assume that an instability of a nonzero spatial 
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wavenumber is for k = 1. We will furthermore co-expand 
the constant input / so as to maintain a fixed value for 
the spatially homogeneous steady state solution R 



1 = 1 



R + eri + 



(14) 
(15) 



where the small parameter e is defined by Eq |13| Addi- 
tionally we define the slow time 



A. Turing Bifurcation: lu = Q, k ^ 



(16) 



The emergence and nature of stationary bumps in rate 
equations have been extensively studied elsewhere, e.g. 
[2H]. We briefly describe this state here for complete- 
ness. The A:*'* spatial Fourier mode of the connectivity is 
given by the critical value Eq. |9] while we assume that all 
other Fourier modes are sufficiently below their critical 
values to avoid additional instabilities. Without loss of 
generality we assume k = 1 here. 

We expand the parameters Ji, / and r as in 
Eqs. |13|14 | l"5l and define the slow time Eq.[l6l The solu- 
tion of Eq|3] linearized about the SU state i? is a spatially 
periodic amplitude which we allow to vary slowly in time, 
i.e. ri = A(T)e^^ + c.c Carrying out a weakly nonlin- 
ear analysis to third order in e leads to the amplitude 
equation 



dTA = r]AJkA + T\A\^A, 



with the coefficients 



1 + D' 



(17) 
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FIG. 3: (Color online) The phase diagram for stationary 
bumps as a function of the zeroth and second spatial Fourier 
modes of the connectivity kernel. The region of bistability 
between the unpatterned and the bump state is shaded. Here 
the critical spatial Fourier coefficient Ji = 3.54. Red fines 
ndicate the boundaries of the SU state (obtained via 



10b I. The functions $ and J{x) as wefl as the input current 
/ and the defay D are taken as in Fig|2] Insets: exampfe 
connectivity patterns corresponding to the vafues of Jo and 
J2 marked by the square and triangfe respectivefy. Note that 
standard Mexican Hat connectivity tends to favor bistability. 



a saddle-node bifurcation for values of Ji below the criti- 
cal value for the Turing instability. Such finite-amplitude 
bumps are therefore bistable with the SU state. In Fig|3] 
the two insets show the connectivity kernel J{x) for pa- 
rameter values given by the placement of the open trian- 
gle (subcritical bump) and the open square (supercritical 
bump) . 

In the phase diagram Figj2] the Turing instability line 
(upper horizontal red line) is shown thin for supercritical, 
and thick for subcritical bumps (here J2 = 0). 



1 + D\1-Jq<^' 2(1 -J2$') 



(18b) 



The nature of the bifurcation (sub- or supercritical) 
clearly depends strongly on the sign and magnitude of 
mean connectivity Jq and the second spatial Fourier 
mode J2. FigjSjshows a phase diagram of the bump state 
at the critical value of Ji = 3.54. The red lines indi- 
cate oscillatory and steady instability boundaries for the 
modes Jq and J2 • Clearly Jo < and J2 < over most of 
the region of allowable values, and the bump is therefore 
supercritical. There is only a narrow region of predomi- 
nantly positive values (shaded region in Figj3]) for which 
the cubic coefficient is positive. This indicates that the 
bifurcating solution branch is unstable. However, neu- 
ronal activity is bounded, which is captured in Eqj3]by a 
saturating transfer function <I>. Thus the instability will 
not grow without bound but rather will saturate, produc- 
ing a finite amplitude bump solution. This stable, large 
amplitude branch and the unstable branch annihilate in 



B. Hopf Bifurcation: a; 7^ 0, fc = 

There is a spatially homogeneous oscillatory instability 
with frequency uj given by Eq |10a[ This occurs for a 
value of the 0*'' spatial Fourier mode of the connectivity 
given by Eq[lObl while we assume that all other Fourier 
modes are sufficiently below their critical values to avoid 
additional instabilities. We expand the parameters Jq, / 
and r as in Eqs. 13|14|15 and define the slow time Eq. 16 
The linear solution has an amplitude which we allow to 
vary slowly in time, i.e. ri = H{T)e^^* + c.c. Carrying 
out a weakly nonlinear analysis to third order in e leads 
to the amplitude equation 



^TH={^l + in)AJoH + {a + il3)\H\^H, 



(19) 



where the coefficients (/i -t- ifi) and (a -I- i(3) are specified 
by the Eqs A4 and |A5| in the Appendix. 

Fig|4] shows a typical bifurcation diagram (in this case 
Ji = 0) for the Hopf bifurcation. Plotted is the ampli- 
tude of the limit cycle as a function of Jq where symbols 
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-0.08 




FIG. 4: The bifurcation diagram for a supercritical Ifopf bi- 
furcation. Shown is the ampHtude of the hmit cycle as a func- 
tion of the 0*'' order spatial Fourier coefficient of the coupling 
J{x). Open circles are from numerical simulation of Eqjsjand 
solid lines show the solution from the amplitude equation, 



Eq 17 The functions $ and J{x) as well as the input current 
/ and the delay D are taken as in Fig[2] 



are from numerical simulation of EqjS] and the lines are 
from the amplitude equation, Eq{T9] 

In the small delay limit {D 0) we can use the asymp- 



totic values Eqs 12 to obtain, to leading order, 



in 



(20a) 



a + ip 



X / (11^-4) ($ )2 



20 



7r<i> 



(ll + 7r) 

10 2 



(20b) 



where we have defined the quantity x = / (8 + 2n^). 

FigjSjshows a comparison of the full expressions for the 
coefficients of the amplitude equation, Eqs A4]|A5 with 
the expressions obtained in the limit D ^ 0, Eqs |20a| 
- |20b| Again, the agreement is quite good, even up to 
D — 1, especially for the real part of the cubic coeffi- 
cient a, which is of primary interest here. The asymp- 
totic expression for the cubic coefficient a, Eq |20b[ in- 
dicates that a subcritical limit cycle should occur for 
<i>"'<i>7(<&")^ > (IItt - 4)/(57r). This provides a sim- 
ple criterion for determining whether or not a particular 
choice of the transfer function can generate oscillations 
which are bistable with the SU state. In fact, it is a 
difficult condition to fulfill given a sigmoidal-like input- 
output function. For example, given a sigmoid of the 



form <I>(x) = a/{l 



), one finds that 



-3/3 a; 



(Ci,")2 



(21) 




FIG. 5: Top: The real part of the linear coefficient ^. Bottom: 
Minus the real part of the cubic coefficient —a. Solid lines are 
from the full expressions Eq |A4)A5l and dotted lines are the 
leading order terms in the small delay limit, Eqs |20a|20b| The 
functions <1> and J{x) as well as the input current / are taken 
as in Fig[2] 



It is straightforward to show that Eqj2Tj is bounded 
above by 1. In fact, -oo < <&'" <&'/(<& ')^ < 1 < 
(11 — 47r) / (Stt) ~ 1.95. Such a nonlinear transfer function 
will therefore always generate supercritical oscillations. 

If the nonlinear transfer function is interpreted as the 
single-cell fl curve, which is common in the literature, 
then we can use the fact that cortical cells operate in 
the fluctuation-driven regime. In particular, the mean 
input current to cortical cells is too low to cause spik- 
ing. Rather, this occurs at very low rates due to fluctu- 
ations in the membrane voltage. Although the fl curve 
for spiking neurons in the supra-threshold regime is con- 
cave down and saturates, in the fluctuation-driven, sub- 
threshold regime the fl curve exhibits a smoothed out 
tail which is concave up. It has been shown that the 
sub-threshold portion of the fl curve of actual cells can 
be well flt by a function of the form ^(x) = Ax'^ , where 
7 > 1, see e.g. K is straightforward to show that 

in this case 



1 



7-r 



which again is bounded between —oo and 1. This again 
rules out subcritical oscillations in the small delay limit. 
Nonetheless, suitable functions <i> for generating subcrit- 
ical oscillations can be contrived, as shown in Fig|6}\. 
Numerical simulation of Eq|3] indeed reveals a subcrit- 
ical bifurcation in this case, see Figj6}3. However, this 
type of transfer function does not seem consistent with 
the interpretation of <i> as a single-cell fl curve, nor with 
that of <i> as a cumulative distribution of activation, i.e. 
a sigmoid. This strongly suggests that delay-driven os- 
cillations in networks of spiking neurons will be generi- 
cally supercritical. This is consistent with the result of a 
weakly nonlinear analysis of a Hopf instability in a net- 
work of integrate-and-flre neurons [4 and with numerical 
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(A) 2 




FIG. 6: A. An example of a function <l>(a;) for which subcrit- 
ical oscillations are possible. The dotted curve indicates the 
range of the function <E>over which oscillations are subcritical. 
B. A bifurcation diagram for subcritical oscillations when the 
function ^{x) is the same as in panel A. Open circles: the 
limit cycle amplitude computed numerically as a function of 
Jo- Here D — 0.1 and the critical coupling is Jo = —15.89. 
The fixed point is held at ii = 0.1 and thus the value of x in 
panel A is close to 0.1 {x + x"^ — 0.1). 



simulations of Hodgkin-Huxlcy neurons ^ . 



C. Turing-Hopf Bifurcation: lj 0, k ^ 

There is a spatially inhomogeneous oscillatory instabil- 
ity with frequency u given by Eq 10a This occurs for a 
value of the fc*'* spatial Fourier mode of the connectivity 



given by Eq |10b| while we assume that all other Fourier 
modes are sufficiently below their critical values to avoid 
additional instabilities. Without loss of generality we as- 
sume k = I. 

We expand the parameters Ji, / and 
Eqs. |13|14|15] 



16 



r as m 
The hn- 



and define the slow time Eq, 
ear solution consists of leftwards and rightwards traveling 
waves with an amplitude which we allow to vary slowly 




FIG. 7: Top: The real part of the cubic coefficient a. Bottom: 
The real part of the cross-coupling coefficient c. Solid lines 
are from the full expressions Eq ]A6|A7l and dotted lines are 
the leading order terms in the small delay limit, Eqs |23a)23b| 
The functions $ and J{x) as well as the input current I are 
taken as in Fig[2] 



rying out a weakly nonlinear analysis to third order in e 
leads to the coupled amplitude equations 

dxA = + in)AJiA +{a + ib)\A\^A +{c+ id)\B\^A, 

(22a) 

drB = (a* - in)AJiB + (a - ib)\B\^B + (c - id)\A\^B, 

(22b) 

where the coefficients (a + ib), (c -I- id) and (fi + ifl) are 
given by the Eqs |A6| |A7| and [A4[ respectively. 

In the small delay limit (D — s- 0) we can use the asymp- 
totic values EqsfT2|to obtain, to leading order. 



a + ib 



x(f + z) / M^")' 

{D^')3 I 1 - $'Jo 2 



c + id 



x(f + ^) / M<i>"r 

(Ll$')3 I l-$'Jo 




(23a) 



(23b) 



Figlr] shows a comparison of the full expressions (solid 
lines) for the real parts of the cubic and cross-coupling 
coefficients a and c with the asymptotic expressions 
above (dotted Hnes). 



Wave solutions and their stability. 



Eqs 22a and 22b admit solutions of the form (A, B) = 
{Ae^^-^ , Be"^^^), where the amplitudes A and B obey 



in time, i.e. ri = A(T)e 



-BiT)e 



-iujt+ix _ 



-C.C.. Car- 



A = iiAJiA+ aA^ + cB^A, 



B = nAJiB + aB^ + cA^B. 



(24a) 



(24b) 
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Traveling waves: Leftward and rightward traveling 



waves in Eqs 24a and 24b are given by {Atw,0) and 
{0,Atw) respectively, where Atw = — /iAJi/a. The 
stability of traveling waves can be determined with the 
ansatz {A,B) = {Atw,0) + {SA,SB)e^*. The resulting 
eigenvalues are Ai — — 2/iAJi and A2 = — /lA Ji(c/a— 1). 



Standing waves: Standing waves in Eqs |24a| and |24b| 
are given by {Asw,-^sw), where Atw = —fJ,AJi/{a+c). 
The stability of standing waves can be determined with 
the ansatz {A,B) = {Asw,Asw) + {SA,SB)e^*^. 
The resulting eigenvalues are Ai = — 2/iAJi and 
A2 = -2/^AJi(a- c)/(a + c). 

The existence and stability of small-amplitude waves 
as described above is completely determined by the val- 
ues of the cubic and cross-coupling coefficients a and c. 
This is illustrated in Fig|8] where the parameter space 
is divided into five sectors. In each sector the type of 
solution which will be observed numerically is indicated 
where known, and otherwise a question mark is placed. 
Illustrative bifurcation diagrams are also given. Specifi- 
cally, in the region labeled 1 (red online), the SW solution 
is supercritical and unstable while the TW solution is su- 
percritical and stable. TW will therefore be observed. In 
the region labeled 2, the SW solution is supercritical and 
unstable while the TW solution is subcritical. Finite- 
amplitude TW are therefore expected to occur past the 
bifurcation point. In the region labeled 3, both solution 
branches are subcritical, indicating that the analysis up 
to cubic order is not sufficient to identify the type of 
wave which will be observed. In the region labeled 4, 
TW are supercritical and unstable while SW are subcrit- 
ical. Finite amplitude SW are therefore expected past 
the bifurcation point. In the region labeled 5, the TW 
solution is supercritical and unstable while the SW so- 
lution is supercritical and stable. SW will therefore be 
observed. 

In the small delay limit, we find, from Eqs 23a|23b 
that 



X 



2 (£1$)= 



(25) 



It is clear from Fig|8] that the nature of the solution 
seen will depend crucially on the sign of the second term 
of the right-hand side of Eq|25) In particular, the di- 
agonal a — c divides the the parameter space into two 
qualitatively different regions. Above this line TWs are 
favored while below it SWs are favored. In the small de- 
lay limit, Eq|25] indicates that the balance between the 
third derivative of the transfer function $ and the sec- 
ond spatial Fourier mode of the connectivity kernel will 
determine whether TW or SW are favored. 

For sigmoidal transfer functions, the third derivative 
changes sign from positive to negative already below the 
inflection point, while for expansive power-law nonlin- 
earities, which fit cortical neuronal responses quite well 
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FIG. 8: The existence and stability of traveling and standing 
waves as a function of the cubic and cross-coupling coefficients 
a and c. In each sector of parameter space a representative 
bifurcation diagram is shown. Supercritical (subcriticaf) so- 
futions are shown growing from feft to right (right to feft). 
Stable (unstable) sofutions are given by sofid (dashed) lines. 
Also indicated in each sector is the type of sofution which 
wiff be seen numericaffy. A question mark is pfaced wherever 
the type of stabfe sofution cannot be determined through a 
weakfy nonlinear anafysis. 



in the ftuctuation-driven regime, the third derivative is 
positive if the power is greater than 2 and is negative 
otherwise. The contribution of this term therefore will 
depend on the details of the neuronal response. In sim- 
ulations of large networks of conductance-based neurons 
in the fluctuation-driven regime in which ,Ji was zero, 
the standing wave state was always observed, indicating 
a > [1, 2 . 

The phase diagram for J2 = 0, Figj2] clearly shows the 
dominance of the SW solution, indicating <i> > for 
the parameter values chosen. Specifically, for values of 
Jo < —6.3, supercritical standing waves are stable, see 
region 5 in Fig|8] Figj9]A and Figj9^ show supercritical 
SW patterns for Jo = —40 and Jo — —9 respectively. 
For —6.3 < Jo < —2.6 TW are supercritical and un- 
stable while SW are subcritical, see region 4 in Fig|8] 
An example of subcritical SW is shown in Figj9p. For 
-2.6 < Jo < 3.58 both SW and TW are subcritical, see 
region 3 in Fig|8] Numerical simulations reveal TW pat- 
terns in this region, see an example in Figj9jl). In the 
region where SW are subcritical there is a small sliver in 
( Jo, Ji) where the SW state is bistable with a TW state 
(TW/SW in the phase diagram). This TW branch most 
likely arises in a secondary bifurcation slightly below the 
Turing-Hopf bifurcation line. There is also a small region 
of bist ability between large amphtude TW and the spa- 
tially uniform high activity state (TW/HA in the phase 
diagram Figj2]). 

Thus for <i> > and with a simple cosine connectivity, 
SW arise for most values of Jo. However, adding a non- 
zero J2 can lead to the TW solution winning out. The 
phase diagram of wave states as a function of Jq and J2 
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FIG. 9: Examples of wave solutions from numerical simulation 
of Eqjs] The functions $ and J{x) as well as the input current 
/ and the delay D are taken as in Fig[2) with Ji = -120. A. 
Supercritical standing waves: Jo = —40 and 5 units of time 
are shown. B. Supercritical standing waves: Jq — —9 and 
5 units of time are shown. C. Subcritical standing waves: 
Jo = —5 and 40 units of time are shown. D. Subcritical 
traveling waves: Jo = and 5 units of time are shown. 



is shown in Fig|lO}\. In FigjTOjV, the hght shaded region 
indicates values of Jq and J2 for which TW are expected, 
whereas SW are expected in the dark shaded region. In 
the unshaded region, both TW and SW are subcritical 
and the solution type is therefore not determined by the 
analysis up to cubic order. These regions, delimited by 
the solid lines, were determined by numerically evaluat- 
ing the real parts of the full expressions for the cubic 
and cross-coupling coefficients, Eqs |A6|X7l Each region 
is furthermore numbered according to the existence and 
stability of the TW and SW solution branches as shown 
in FigjS] The dashed lines show the approximation to the 
solid lines given by the asymptotic formulas Eqs |23a||23b[ 
The set of allowable values for Jq and J2 is bounded by 
the conditions for a steady or oscillatory linear instability 
Eqs |9p0bl These stability condit ions are shown by the 
horizontal and vertical bounding lines (red online). All 
parameter values are as in Fig|2] 

From Fig |10| we can now understand the discrepancy 
between the analytical results in [T] using a rate equation 
with a linear threshold transfer function, which predicted 
TW, and network simulations, which showed SW. Specif- 
ically, given a nonlinear transfer function with 4> > 0, 
then with a simple cosine coupling SW are predicted over 
almost the entire range of allowable Jq (dark shaded re- 
gion for J2 — 0). The nonlinear transformation of inputs 
into outputs is thus crucial in determining the type of 
wave solution. The choice of a threshold linear trans- 
fer function results in the second and all higher order 
derivatives being zero. In this sense it produces degener- 
ate behavior at a bifurcation point, and by continuation 
of the solution branches, in a finite region of the phase 
diagram. 
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FIG. 10: (Color online) A. Phase diagram for waves as a func- 
tion of the zeroth and second spatial Fourier coefficients of the 
connectivity kernel. The dark-shaded region indicates SW, 
whereas the light shaded region indicates TW. Red lines indi- 
cate boundaries for primary instabilities of Jo and J2 given by 
Eqs |9|10"bl Solid stability lines for waves are from Eqs |A6|A7| 
whi le the da shed line are from the asymptotic expressions, 
Eqs |23a|23bl Here Ji = -58.4. The function $ as well as 
the input current / and the delay D are taken as in Fig[2] 
Insets: example connectivity patterns corresponding to the 
values of Jo and J2 marked by the square and triangle re- 
spectively. B. The same phase diagram as in A, now showing 
where various types of 'difference-of-Gaussian' connectivities, 
Eq. |26| would lie. Each dotted line indicates the border of a 
region in which the standard deviations of the excitatory and 
inhibitory connectivities are below a certain threshold (0.7, 
1.0 and 1.5 respectively, see text). Relatively narrow connec- 
tivities compared to the system size will always generate TW 
solutions. See text for details. 



We have shown that varying Jq can change the nature 
of the bifurcation, e.g. supercritical to subcritical, while 
varying J2 can switch the solution type, e.g. from SW to 
TW. As an example of a functional form of connectivity 
motivated by anatomical findings, e.g. |31j . we consider 
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a difference of Gaussians, written as 

J(x) = — i^e"^ - — i^e"^. (26) 
In this case, one finds that the Fourier coefficients are 
Jfc - Jee-^-'-^'/2/(fc,ae) - J,e-'=''^'/V(fc,a.)- (27) 

where /(fc,CTe,,) = Re[Erf((7r/cre,, + iA:^)/%/2)]/7r. Once 
Ji has been fixed at the critical value for the onset of 
waves, from Eq 27 it is straightforward to show that 
Jo = — pJ2 + 1 where both p and q are constants which 
depend on and cr^, the widths of the excitatory 
and inhibitory axonal projections respectively. Thus 
a difference-of-Gaussian connectivity, constrains the 
possible values of Jq and J2 to lie along a straight 
line for fixed connectivity widths. This is illustrated 
in Fig jlOp where three dashed lines are superimposed 
on the phase diagram, corresponding to the values 
de,, = (1.5,1.49); ^Te,^ = (1,0.99); and a,^, = (0.7,0.69). 
Each of these lines is bounding a region to the left where 
CTe and <Ji are less than 0.7, 1.0 and 1.5 respectively. 
Given periodic boundary conditions with a system size 
of 27r, a Gaussian with a — 1.5 is already significantly 
larger than zero for x = tt or — tt. Thus, restricting 
ourselves to Gaussians which essentially decay to zero at 
the boundaries means that TW will always be observed. 
The same holds true for qualitatively similar types of 
connectivity. 



2. Waves in Network Simulations 

Our analytical results concerning waves from the rate 
equation EqjS] predict that a connectivity with a suffi- 
ciently strong second Fourier component with a negative 
amplitude will lead to traveling waves. We can test this 
prediction in simulations of networks of spiking neurons. 
In doing so we use the same general network as in ear- 
lier work [Tl, 15] but with modified connectivity. Briefly, we 
consider two populations of iV = 2000 conductance based 
neurons each, arranged on a ring. We use a Hodgkin- 
Huxley formalism with one somatic compartment. Ac- 
tion potentials are shaped by Na and K. For details of 
the model see [52] . The probability of connection from a 
neuron in population A{= E, I) to a neuron in population 
B{= E, I) is pb A — Po^ + Pi^ cos r -I- ^ cos 2r, where 
r is the distance between the two neurons. Here we take 

P. 



2000^ 



BA 



pf, i.e. the excitatory and inhibitory input is the 
same for both excitatory and inhibitory neurons on av- 
erage. This leads to one effective population of neurons, 
the appropriate condition for comparing to EqjSj Synap- 
tic currents are modeled as Isyn.A — ~9As{t)(V — Va), 
A G E,I where s{t) is modeled as the difference of ex- 
ponentials with time constants equal to one and three 
milliseconds respectively, V is the post-synaptic voltage 
and Va is the reversal potential and is mV and —80 
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FIG. 11: Raster plots of spiking activity in a network of 
conductance-based neurons. A. Standing waves. Parame- 
ters are the same as in Fig. 3 [1], right column, third down. 
Pe = 0.2, Pi = 0.2 + 0.2 cos r. B. Traveling waves, pe = 0.2, 
pi = 0.2 0.2 COST- -I- 0.1 cos 2r. gs = 0.01, gi = 0.028, 
ge:.t = 0.001, u^^t = 5000. 



mV for excitatory and inhibitory synapses respectively. 
A post-synaptic response s{t) is initiated whenever a pre- 
synaptic action-potential reaches mV, with no delay. 
All the neurons are subjected to an external excitatory 
input comprised of synaptic activations {Isyn,E as defined 
above) with Poisson statistics at a rate i/gajt- Fig|TT]shows 
the results of the simulations. As predicted the addition 
of the second spatial Fourier component to the inhibitory 
connections converts SW into TW. 



IV. BIFURCATIONS OF CODIMENSION 2 

For certain connectivity kernels we may be in the vicin- 
ity of two distinct instabilities. This is the case for certain 
Mexican hat connectivities (OU and SB) and certain in- 
verted Mexican hat connectivities (OU and SW/TW). 
Although two instabilities will co-occur only at a single 
point in the phase diagram Fig|2] i.e. Jq and Ji are 
both at their critical values, the competition between 
these instabilities may lead to solutions which persist 
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over a broad range of connectivities. This is the case 
here. We can investigate this competition once again us- 
ing a weakly nonhnear approach. The main results of 
this section are summarized in table HI 



A. Hopf and Turing-Hopf bifurcations 

Here we consider the co-occurrence a spatially homo- 
geneous oscillation and a spatially inhomogeneous oscil- 
lation (OU and SW/TW), both with frequency u) given 
by Eq |10a[ This instability occurs when the zeroth and 
fc*'* spatial Fourier mode of the connectivity both sat- 
isfy the relation, Eq |10b| while we assume that all other 
Fourier modes are sufficiently below their critical values 
to avoid additional instabilities. Without loss of general- 
ity we take fc = 1 for the SW/TW state. 

We expand the parameters Jq, Ji, I and r as in 
13|14|15| 



Eqs. 
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and define the slow time Eq. 
linear solution consists of homogeneous, global oscil- 
lations, leftwards and rightwards traveling waves with 
amplitudes which we allow to vary slowly in time, i.e. 



H{T)e"^* + A{T) 



-,iLjt-{-ix ^ 



B(T)e-'"*+"-Hc.c.. Car- 



rying out a weakly nonlinear analysis to third order in e 
leads to the coupled amplitude equations 



+2{a + if3)[{ 



(28a) 

\A\^ + \B\^)H + HAB], 



drA = {fi + in)AJiA+{a + ib)\AfA+{c + id)\BfA 
+ia + iP)[2\H\^A + H^Bl (28b) 



Oscillatory Uniform (OU): The oscillatory uniform solu- 
tion has the form {H, A, B) = (He'"*, 0, 0) where 



n = 



Lu = [ft - —jAJo- 

The OU state undergoes a steady instability along the 
line 



AJi = AJo- 



(29) 



This stability line agrees very well with the results of nu- 
merical simulations of EqJS} see the phase diagram Fig|2] 

Traveling Waves (TW): The traveling wave solution has 
the form [H, A, B) = (0, Atwg'-'^^Q) or (0, 0, Arwe-^'^^), 
where 



A 



TW 



a 

c ^ (^^-^)aji. 

The TW state undergoes an oscillatory instability 
along the line 



AJi ^ —AJo, 



(30) 



with a frequency 

^.(^.(1- A) + (5 -2^) A) A Jo. 

Standing Waves (SW): The standing wave solution has 
the form {H,A,B) = {0, Aswe''^^ Aswe''"^^), where 



drB = {n-in)AJiB + {a-ib)\B\^B + {c-id)\A\^B 
+ia-iP)[2\H\^B + H^A], (28c) 



whe re a + i(3, a + ib and c + id are given by Eqs |A5| |A6| 
and A7 respectively. The overbar in H represents the 
complex conjugate. 



Solution types and tfieir stability 



Eqs |28a|28c| admit several types of steady state 
solutions including oscillatory uniform solutions (OU), 
traveling waves (TW), standing waves (SW) and mixed 
mode oscillations/standing waves (OU-SW). The sta- 
bility of these solutions depends on the values of the 
coefficients in Eqs |28a|28c[ In addition, non-stationary 
solutions are also possible. Here we describe briefly some 
stationary solutions and their stability. For details see 
Appendix 



As 



w 



^jiAJi 
/(« + c) ' 

'n-^^n)AA. 



{a + c) 

An oscillatory instability occurs along the line 



4a 



with a frequency 



(31) 
(32) 

(33) 
AJo- 



a + c b + d-AI3 + 

A stationary instability occurs along the line 

AJi = ^-AJo, (34) 

where 



-k2 + ^Jk\ -Akiks 



2fci 
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il — /i 



(a + c) 



2(120^^^4/32) 
^ (a + c)2 '■ 



[a + c) [a + c) 



k3 = n^ + n\ 

For Eqj3]with the parameters used to generate the phase 
diagram, Figj2]we find that the stationary instabihty pre- 
cedes the oscillatory one and that "if ^ 0.6. This agrees 
well with the numerically determined stability line near 
the co-dimension 2 point, see Figj2j 

Mixed Mode: We can study the mixed mode solutions 
m Eqs |28ap8c| by assuming an ansatz 



(36) 



which leads to four coupled equations, see Appendix [A] 
We do not study the stability of mixed mode solutions in 
this work. 



2. A simple example 

We now turn to a simple example in order to illustrate 
the two main types bifurcation scenarios that can arise 
when small amplitude waves and oscillations interact in 
harmonic resonance. 

i. Bistability: Here we take the parameters 01]. Given 
these parameter values one finds, from the analysis above, 
that the oscillatory uniform state has an amplitude Ti. = 1 
and destabilizes along the line AJk — —1. The standing 
wave solution (traveling waves are unstable, see Figjs]) 
has an amplitude Asw = V~^Jk which undergoes a 
steady bifurcation to the oscillatory uniform state at 
AJk — —1/2. Both solutions are therefore stable in the 
region —1 < AJ^, < —1/2. This analysis is borne out by 
numerical simulation of Eqs |28aj|28c| see Fig 12 1. Solid 
and dotted lines are the analytical expressions for the 
stable and unstable solution branches respectively (red 
is OU and black is SW) . Circles are from numerical sim- 



ulation of the amplitude equations Eqs 28a - 28c 



Note that this scenario is the relevant one for the 
phase diagram shown in Fig|2] That is, we find there is a 
region of bistability between the OU and SW solutions, 
bounded between two lines with slope ~ 0.6 and 1 
respectively. 

a. Mixed Mode: Here we consider the parameters [55] . 
Given these parameter values one finds that the oscilla- 
tory uniform state has an amplitude Ti. — I and desta- 
bilizes along the line AJ^ = —1. The standing waves 
solution has an amplitude Asw — ■\/— A Jfc/8 (traveling 



waves are again unstable) which undergoes an oscilla- 
tory instability at AJ^ = —2. The mixed-mode solution 
is given by 



n 

Asw 



8 + 2A Jfc(2 - cos (/) - sin ( 
8-2(2-cos(/)-sin(?!))2 
A Jfe -|- (2 — cos — sin </)) 
8 - 2(2 - cos 0- sin 



A J/j(l — 4cos — 2 sin < 
—4 cos 0-1- 8 cos 4? . 



2 sin 6 cos < 



(37a) 
(37b) 

I 

(37c) 



It is easy to show that for AJfc —1 the mixed mode 
ampUtudes approach (Ji.Asw) — (1,0) and the phase 
— > 0. The mixed-mode solution thus bifurcates contin- 
uously from the oscillatory pure mode. FigfT2|3 shows the 
corresponding bifurcation diagram where solid and dot- 
ted lines are the analytical expressions for the solution 
branches and symbols are from numerical simulation of 
Eqs |28a| - |28c| As AJ^ increase from the left we see that 
the SW solution indeed undergoes an oscillatory insta- 
bility at AJfc = —2 leading to an oscillatory mixed- mode 
solution indicated by small circles (the maximum and 
minimum amplitude achieved on each cycle is shown). 
This oscillatory solution disappears in a saddle-node bi- 
furcation, giving rise to a steady mixed-mode solution 
whose amplitude is given by Eq 37a- 37c This steady 



mixed-mode solution bifurcates from the pure oscillatory 
mode at AJ^ = — 1 as predicted. 



3. Summary 

The interaction between the oscillatory uniform state 
and waves may lead to mixed mode solutions or bista- 
bility. The OU state always destabilizes along the line 
Ji — Jo, irrespective of parameter values or the choice of 
$ or J{x). This result from the weakly nonlinear anal- 
ysis, agrees with numerical simulations of Eqj3] over the 
entire range of values of Jq and Ji used in the phase dia- 
gram, Figj2] and appears to be exact. Depending on the 
value of J2 , supercritical TW or supercritical SW will be 
stable near the codimension 2 point. In the case of TW, 
the slope of the stability line is one half the ratio of the 
cubic coefficient of waves to that of oscillations. In the 
small delay limit this expression can be simplified to 



a 

2a 



4 ( (ll7r^4) ($")^ 
20 *' 



4 



(38) 



which depends only on shape of the transfer function 
For the parameter values used in the phase diagram Figj2] 
this yields a line with slope close to one half. Thus TW 
and OU are expected to be bistable in the wedge between 
AJi = AJo/2 and AJi = AJq- In the case of SW, the 
slope of the stability line is a complicated function of the 
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FIG. 12: Two typical bifurcation diagrams for the case 
of harmonic resonance between small-amphtude oscillations 
and small-amplitude standing waves. A: Here oscillations 
and standing waves are bistable for —1 < AJi < —1/2. 
AJo = -1, a = -1, a = -1, 6 = c = d = /3 = n = 0, 
= — 1. B: Here the standing wave solution loses stabil- 
ity to an oscillatory mixed-mode solution at AJi = —2. At 
AJi ~ —1.75 a steady mixed-mode solution arises in a saddle- 
node bifurcation and continuously approaches the oscillatory 
pure- mode solution at AJ^ = —1. Parameters are a = — 8, 
/3 = 1, Q = —1, PL = —1, h = c = d = VL = Q. The phase (/) of 
the mixed-mode solution is not shown. 



shape of <i> and the second Fourier coefficient J 2 ■ For the 
parameter values used in the phase diagram Figj2] the 
slope is close to 0.6. Therefore the OU and SW states 
are bistable in the wedge between AJi = O.6AJ0 and 
AJi = AJq. 



B. Hopf and Turing bifurcations 

We consider the co-occurrence of two instabilities: a 
spatially homogeneous oscillation and a spatially inho- 
mogeneous steady solution. This occurs when the zeroth 
spatial Fourier mode of the connectivity satisfies the re- 
lation, Eq |10b| and the kth spatial Fourier mode satisfies 
Jfc = 1/$ , while we assume that all other Fourier modes 
are sufficiently below their critical values to avoid addi- 
tional instabilities. Without loss of generality we take 
fc = 1 for the Turing instability. 

We expand the parameters Jq, Ji, / and r as in 



Eqs. |13|14|15| and define the slow time Eq. (Te] The 
linear solution consists of homogeneous, global oscilla- 
tions and stationary, spatially periodic bumps with am- 
plitudes which we allow to vary slowly in time, i.e. 
n = i7(r)e*'^* + A(r)e*'=^ + c.c. Carrying out a weakly 
nonlinear analysis to third order in e leads to the coupled 
amplitude equations 

drH = {^i + m)/^jQH + {a + il3)\H\'^H 

+{n + iK)\A\^H, (39a) 

dTA^f]A.hA + V\A\^A + (j\H\'^A, (39b) 
where /i + a + ij3, f], T , k + iK and a are given by 



Eqs A4 A5 18a 18b A15 and A16 respectively. 



1, Solution types and their stability 



Steady state solutions to Eqs 39a and 39b include pure 
mode OU, pure mode SB and mixed mode solutions (OU- 
SB). We look at the stability of the OU and SB solutions 
in turn for the general case and then look specifically 
at the case of small delay in Eq|3] Since the coupling 
in Eqs |39a| - |39b| is only through the amplitudes we can 
simphfy the equations by taking {H,A) — {He"^^ , Ae^'^) 
which yields 



A = T]A,hA + TA^ + crH^A. 



(40a) 



(40b) 



Uniform oscillations (OU): The oscillatory uniform solu- 
tion has the form (7Y, 0) where 



mAJq 
a 



The linear stability of this solution can be calculated with 
the ansatz 

{n,A) ^ {n + Sne^\SAe^'), 
which yields the two eigenvalues 

Xh = -2/iAJo, 



\a = ?7A Ji + —AJo. 

a 

If we assume a supercritical uniform oscillatory state then 
the first eigenvalue is always negative, while the second 
becomes positive along the line 



AJi = ^AJo, 
indicating the growth of a bump solution. 



(42) 
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For EqjS] with the parameters used to generate the 
phase diagram Figj2] we find from Eq 42 that the OU 
state destabilizes along the line AJi ^ — 0.026A Jq- 

Stationary Bump (SB): The stationary bump solution 
has the form (0,.4) where 



A 



ttAJi 

r 



The linear stability of this solution can be calculated with 
the ansatz 

which yields the two eigenvalues 

\a = -277AJ1. 

If we assume a supercritical stationary bump state then 
the second eigenvalue is always negative, while the first 
becomes positive along the line 



AJi = ^AJo, 



(44) 



indicating the growth of uniform oscillations. 

For EqjS] with the parameters used to generate the 
phase diagram Fig|2] we find from Eq 44 that the SB 
state destabilizes along the line AJi ~ — 0.144AJo. 

Mixed Mode (OU-SB): The mixed-mode solution sat- 
isfies the following matrix equation 



which yields 



a K 

G r 



re 

A^ 



\ '?AJi 



A^ = 



-fiTAJo + tjkAJi 
aT — an 



IIuAJq — rjaAJi 
aT — UK 



We do not study the stability of the mixed-mode solution 
here. 



2. A simple example 

We once again illustrate the scenarios of bistability 
and mixed-mode solutions with a simple example. 

i. Bistability: Here we take the parameters |36j . In 
this case, the limit cycle has an amplitude 7i = 1 and 
undergoes an instability at A Ji = 2. The bump solution 
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FIG. 13: Typical bifurcation diagrams for the competition 
between bumps and global oscillations, [s, — —1, A Jo = —1, 
a = -1, = 1, r = -1. A: K = -2, a = -2. The limit cycle 
and bump solutions are bistable in the range 1/2 < AJi < 2. 
B: A mixed-mode solution is stable in the range 1/2 < AJi < 
2. K = —0.5, a = —0.5. Symbols are from simulation of the 
amplitude equations Eqs |40a|40b| while lines are the analytical 
expressions. 



has an amplitude A = ^/AJi and becomes unstable at 
AJi = 1/2. The oscillatory and bump solutions are 
therefore bistable in the range 1/2 < AJi < 2. This is 
borne out in numerical simulations of Eqs |40a||40b| see 
Fig 13 Symbols are from numerical simulation (cir- 



clesdimit cycle, squares :bump), while lines arc analytical 
solutions. 



ii. Mixed-mode: Here we consider the parameters |37j . 
In this case, the limit cycle has an amplitude Ti = 1 
and undergoes an instability at AJi = 1/2. The bump 
solution has an amplitude A = ^/AJi and becomes un- 
stable at AJi = 2. The mixed-mode solution is stable in 
the range 1/2 < AJi < 2 and has amphtudes TLmm — 
2^(1- AJi/2)/3 and A = 2^(AJi - 1/2)3. The corre- 
sponding bifurcation diagram is shown in Fig |13p where 
symbols are from simulation of Eqs |40a|40b| and lines are 
the analytical results. 
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3. Summary 

The interaction between the SB and OU states can 
lead to one of two scenarios. Either there is a region 
of bistability between bumps and oscillations, or there 
is a mixed-mode solution which, near the codimension 
2 point at least, will consist of bumps whose amplitude 
oscillates in time, i.e. oscillating bumps (OB). 

In the limit of small D the instability lines for the OU 
and SB states in the vicinity of the codimension 2 point 
are given by the equations 



AJi 



AJi 



D 



(ll7r-4) {<!> y 
20 $' 



(46) 



-D 



V 2 2(1- Ja*')^ 



AJo. (47) 



respectively. The slope of both of the stability lines is 
proportional to D, indicating that in the small D limit 
any region of bistability or mixed mode solution will be 
limited to a narrow wedge close to the Jq axis. Which 
scenario will be observed (bistability or mixed-mode) de- 
pends on the particular choice of <i> and the value of the 
second spatial Fourier mode of the connectivity J2. For 
the parameters values used to generate the phase dia- 
gram Figj2] the slopes are ^ —0.026 and ^ —0.144 for 
the OU and SB stability lines respectively, indicating a 
mixed mode solution. 



V. CONCLUSIONS 

In this paper we have studied the nature and stability 
of solutions arising via a primary instability of a fixed 
point in a rate equation with fixed delay, EqjSj The four 
possible primary instabilities are: i. - a steady rate in- 
stability for sufficiently strong recurrent excitation which 
leads to a high activity state, ii. - a Turing instability 
for Mexican hat connectivity which leads to stationary 
bumps, iii. - a Hopf instability for sufficiently strong re- 
current inhibition which leads to an oscillatory uniform 
state, iv. - a Turing-Hopf instability for strong local in- 
hibition and longer-range excitation (inverted Mexican 
hat) which leads to waves. We have focused on the oscil- 
latory instabilities which arise only for a non-zero delay, 
i.e. they are delay-driven. The instability mechanism 
due to the delay can be understood intuitively. We need 
only assume that the coupling in the system is predomi- 
nantly inhibitory and that there is a fixed delay D in the 
interactions between neurons. In this case, if the state of 
the system is perturbed at a time t, causing an increase 
in activity, then this increase in activity will generate a 
corresponding decrease in activity after a time ~ D. This 
decrease in activity leads to an increase in activity again 
after a time D. In this way oscillations can emerge. 



This argument does not tell us how strong the inhibition 
must be to maintain oscillations, and this will depend on 
the details of the system. Additionally it is clear that this 
mechanism is only valid for inhibitory and not excitatory 
feedback. 

An equation of the form Eqj3] with a threshold linear 
transfer function and cosine connectivity was studied al- 
ready in [ij [2] . The simple choice of transfer function and 
connectivity allowed the authors to construct a phase di- 
agram of dynamical states to a large extent analytically. 
Many of the dynamical states predicted by this analysis 
were subsequently confirmed through simulations of re- 
currently coupled spiking neurons [T, ^ . One discrepancy 
concerned the Turing-Hopf instability which in the rate 
equation lead to traveling waves (TW) and in the net- 
work simulations to standing waves (SW). Here we have 
shown that including a transfer function with an expan- 
sive nonlinearity in the rate equation correctly predicts 
a SW state given a simple cosine connectivity. We addi- 
tionally predicted that a TW state should be seen once 
a sufficiently large and negative second spatial Fourier 
mode was included in the connectivity. This was borne 
out in network simulations of Hodgkin-Huxley neurons, 
see Fig |ll[ This explains the discrepancy in the previ- 
ous work but, in fact, makes the much more general point 
that the specific form of the transfer function (or fl curve) 
and the synaptic connectivity may be important in deter- 
mining the dynamical state observed during spontaneous 
activity in neuronal networks. We have shown that the 
same principle holds true regarding the nature of dynam- 
ical states near a bifurcation point, e.g. whether or not 
the solution is subcritical or supercritical or which type 
of solution will be observed if several exist and must com- 
pete as is the case with SW and TW states. 

In general the nature of a dynamical state in EqjSjnear 
one of the primary instability lines depends on a compli- 
cated combination of the first three derivatives of the 
transfer function at the fixed point and the first three 
spatial Fourier coefficients of the synaptic connectivity. 
We have derived these expressions and their asymptotic 
approximations in the limit of small delay. Despite the 
complexity of these expressions we have tried to answer 
two specific qualitative questions: 1 - are delay-driven 
oscillations in general supercritical or subcritical?, 2 - 
given realistic patterns of synaptic connectivity which is 
the most likely wave state, TW or SW? The nature of the 
oscillatory uniform state depends on the balance of the 
first three derivatives of the transfer function. We find 
that whether the transfer function represents the fraction 
of active neurons |6j or the fl curve of a cortical cell driven 
by noisy inputs [23 EU] , the oscillations are likely to be 
supercritical. This agrees with numerical simulations of 
Hodgkin-Huxley neurons (not shown here) as well as with 
the amplitude equation derived analytically for a network 
of integrate-and-fire neurons [4]. We also find that for 
patterns of synaptic connectivity which decay in space, 
e.g. Gaussian, the TW state is favored for small delay. 
The reason for this is that the competition between SW 
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Bifurcations of Codimcnsion 2 



Solution types calculated 



Instability boundaries 



Hopf and Turing-Hopf 



Hopf and Turing 



Oscillatory Uniform 
Travelling Waves 
Standing Waves (oscillatory) 
Standing Waves (stationary) 
Mixed-Mode 
Oscillatory Uniform 
Stationary Bump 
Mixed-Mode (OU-SB) 



AJi = AJo 

AJi = a/ (2a) AJo 

AJi = (a-Hc)/(4a)AJo 

AJi = *AJo 

Not calculated 

AJl = (^Cr)/(77Q)AJo 

AJi = {fir)/{7jn)AJo 
Not calculated 



TABLE I: Some existing d ynamic al s tates that are present close to the codimension 2 bifurcations, and their corresponding 
instability boundaries Eqs. 29 30 33 34 42 44 (except for the Mixed Mode solutions). 



and TW in the small delay limit depends strongly on the 
sign and strength of J2 , the second spatial Fourier mode 
of the connectivity, see Eq 25 At the Turing-Hopf bifur- 
cation point, Ji, the bifurcation parameter, is equal to its 
critical value which is large and negative. The value of J2 
is constrained by Ji to relatively large and negative val- 
ues as well as long as the connectivity is narrow enough 
compared to the width of the system, i.e. it should go 
to zero at the boundaries, in this case at —tt and tt, see 
FigjlO] Large and negative values of J2 favor the TW 
state by Eq|25) Both of these findings are dependent on 
the parameter D in EqjS] being small. As argued in the 
introduction this seems reasonable since D represents the 
ratio of the fixed delay in interactions to the membrane 
time constant which are on the order of 1ms and 10ms 
respectively. 

Finally, we have tried to emphasize the importance 
of fixed delays in shaping the dynamics described by 
Eq|3] and by extension in networks of spiking neurons. 
Nonetheless both fixed and conduction delays are present 
in neuronal systems and are roughly of the same order 
of magnitude in a small patch of cortex of ~ 1mm in ex- 
tent. It remains to be studied how these delays interact 
to shape patterns of spontaneous activity. 
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APPENDIX A: AMPLITUDE EQUATIONS 

In this appendix we outline the calculation of the am- 
plitude equations which describe the slow temporal evo- 
lution of the various instabilities near their respective 
bifurcations. 



1. General framework for the weakly nonlinear 
calculation: Codimension 1 bifurcations 

Here we briefly describe the general framework for the 
weakly nonlinear calculation for the Turing, Hopf and 
Turing-Hopf bifurcations. We use a multiple-scales ap- 
proach which takes advantage of the fact there is a near- 
zero eigenvalue in the vicinity of a bifurcation which is 
responsible for the slow temporal evolution of the critical 
eigenmode. This is a standard method, see e.g. |33j . 

For simplicity we first rewrite Eq|3]as 



-r-K$((Jr) +/), 



(Al) 



where (fg) = ^ J^^ dyf{y- x)g{y,t- D) . We study the 

stability of the steady state solution R = ^(^JqR + , 

where J{x) — Jq + 2 Xl^i cos nx. We expand the 
rates, the connectivity and the input current as 

r{x,t) — R + eri{x,t,T) + er2{x,t,T) + . . . , 
J{x) = J{x) + e^l^J{x), 
I = J-He^AJ, 

where the small parameter e is defined by the distance 
from the critical value of the connectivity. Plugging these 
expansions into Eq XT yields 



(£ + e2/:2)(eri + e^ra + ...) = e'A^2(ri) + e^TVaC'^i, ^2), 
where 

Cr — dtr + r—{Jr), 
C2r = driJr) ~ (AJr), 

n 

N2 - ^{Jri)\ 

N3 = $"(Jri)(Jr2)-f ^(Jri)3. 

We now collect terms by order in e. At first order we 
have 

^{e) : Cri = 0. 
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This equation gives the Hnear dispersion relation Eq|8] 
The values of the connectivity and input current for 
which it is satisfied are J{x) — J{x) and 1 = 1. At 
second order we obtain 

79(e2): Cr2=N2{n). 

The second order solution r2 is the particular solution 
of this linear difl^erential equation. And finally, at third 
order 

d{e^): Cr3 = N3{n,r2) - C2ri. 

At this order secular terms arise which have the same 
temporal and/or spatial frequency as the linear solution. 
In order for the above equation to have a solution, these 
terms must therefore be eliminated, yielding the desired 
amplitude equation for the instability. 



a. Steady Bifurcation: uj = 0, k = 

For completeness we include here the derivation of the 
amplitude equation for the transcritical bifurcation. 

The 0*'* spatial Fourier mode of the connectivity is 
given by the critical value Jq = ^ , while we assume that 
all other Fourier modes are sufficiently below their critical 
values to avoid additional instabilities. We expand 



Jo — Jo 
I = U 



f eA Jo, 



(A2) 



= i? + eri + e^ra 



where the small parameter e is defined by Eq |A2| We 
define the slow time T = et. The linear solution is a 
spatially homogeneous amplitude which we allow to vary 
slowly in time, i.e. ri = -^iT)- Carrying out a weakly 
nonlinear analysis to second order in e leads to the normal 
form for a transcritical bifurcation given by 

drA = ?]AJoA + jA^, 



1] 

1 



1 + £)' 
2(1 + D) 



(A3) 



o Numerics 

— Newton-Raphson 

— Amplitude eqs. 



'0 



FIG. 14; (Color online) Bifurcation diagram for the steady 
instability. Open circles: numerical simulation of Eq[3] Red 
Lines: amplitude equation solution from Eq |A3[ Black lines: 
steady-state solution of Eq|3] using a Newton-Raphson solver. 
Solid lines indicate stable solutions and dotted lines unstable 
ones. $(a;) = x^g'^gx where a = 1.5 and /3 = 3. J(x) = 
Jo + Ji cosx where Ji = Q. The input current / is varied so 
as to keep the uniform stationary solution fixed at _R = 0.1. 



solution are 

N2 

r2 



j2( A2^2ikx 

r22e"""- -f c.c. + r2o 



$ J|(A2g2^fcx^^^^^^2|yl|2)/2, 

lihx 



^"Jl 



2(1- J2fe$')' 

1 - Jo*' ■ 



i?(e^): The nonlinear forcing at cubic order is 

iVg = JkJoAr20 + ^" JkJ2kAr22 

+ Jl\A\'^A/2)e'^'' + C.C. + ... 

Eliminating all terms of periodicity e*''^ at this order 
yields the amplitude equation, EqjlT] 

dTA = riAJuA + T\A\''A, 

with the coefficients 



■n 



73 / 



J2k{^") 



ll\2 



b. Turing bifurcation 



l + D\l-Jo$' 2(1-J2fe$') 



'd{e): The solution to the linear equation is spatially pe- 
riodic with slowly varying amplitude A, 

ri = A{T)e'^'' + c.c. 



c. Hopf bifurcation 

t?(e): The solution to the linear equation is a time peri- 
odic function with slowly varying amplitude H 



-d^e^): The nonlinear forcing and resulting second order 



n = iJ(T)e^"* + c.c. 
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i?(e^): The nonlinear forcing and resulting second order 
solution are 

iV2 = $"jo2(i?2e2»^(t-£') + + 2|iJp)/2, 



r2 = r22e2*'^* + c.c. + r2o, 



" 72 



$ J, 



r-22 



r2Q 



2(2iw + 1 - $'Joe-2»'^-D) 



1 - Jo$' 

'i?(e^): The nonlinear forcing at cubic order is 



d. Turing-Hopf bifurcation 



19(6): The solution to the linear equation are two sets 
of periodic waves, one left-traveling with slowly varying 
amplitude A and the other right-traveling with slowly 
varying amplitude B 



+ Jo^|Apv4/2)e"^(*-°) + c.c. + . . 



Eliminating all terms of periodicity e at this order _ A(r)e'^ -I- _B(r)e'^ -I- c c 

yields the amplitude equation, Eq |19[ 

DtH^ {^i + in)AJaH + {a + if3)\H\^H, 

with the coefficients 

1 -t- 1^(1 + lUj) 
a + i(3= — — X (A5) 

f J^(^ ")'^ J^{^ ")'^e~'^^'^^ J^^ "\ The nonhnear forcing and resulting second order 

°^ /- + — ^ '—^ + -2 . solution are 

ll-$'Jo 2(2icj+ 1 -$'Joe-2«^-0) 2 j 



+2(|Ap-H|Bp))/2, 
r2 = r2^e2'^ + r^^e'^+'> + r^^e'^+^ + r2^e2^-f c.c. + r2o, 

2(2ic^ + l-$'J2fee-2-0) 
r,^ - ^"-^fc AB 



2 



72 

2(-2icji? + l- J2fc$'e2»-0) 



?'20 



<i>"j2(|A|2+|B|2) 



1 - Jo$' 

'i9(e"^): The nonlinear forcing at cubic order is 

N3 = JfcJoAr2o + JfeJ2feAr2v, + <i>"jkJoBr^^ + ^" JkJ2kBr^^ + $"'|ApA/2 + \B\^ A)e*-''^'' 

' lpc/> 



JfeJo^r^^ + $ JfcJ2feAr^0 + JfcJoBrao + $ JkJ2kBr2^ + |S|2b/2 + $ " |A|2s)e'^+''^^ 



Eliminating all terms with dependencies e*^ and e*"^ yields the two coupled amplitude equations Eqs 22a and 22b 



dxA = {n + iQ)\JkA+{a + ih)\A\'^A+{c+id)\B\'^A, 
OtB = {n-in)AJkB+{a-ib)\B\'^B + {c-id)\AfB, 



with the coefficients 



l + WJfee-''^^ \^l-$'Jo 2(2iw + l-$'J2fee-2''^^) 2 J' 
1 + D^'Jke-i'^'^ \ 1 - Jo 2ia; + 1 - Joe-^*'^^ 1 - Jsfe 



2. Codimension 2 bifurcations 

a. Double zero eigenvalue: Hopf, Turing-Hopf 

■!?(e): The solution to the Unear equation are periodic oscillations and traveling waves 

n = H{T) + e*'^* + A(T)e*'^*+''^^ + b{T) + c.c. = iJ(T)e*'^* + A{T)e''' + B{T)e't' + c 

i?(e^): The nonlinear forcing and resulting second order solution are 

+2^Se'^+'^-2»a.i3 ^ g2^24>+2tu,D _^ _^ 2(| + |B|2))/2 

+^"joJk{HAe^"^^*-^'>+" + i5-5e-2*"(*-D)+^^ + ij-^e" + iJBe'^ + c.c), 



' HB 



+rHAe" + THBe''-' + c.c. + rao, 

^ JoJk -2iuD i 

2iw + l-$'jfce-2i'^-c> 
-2zw + l-$'Jfce2»'^-° 



$ JoJfe fV . 

thb = :;fj^HB, 

1 - 4> Jfc 



r20 



1 - Jo$' 

t?(e^): The nonlinear forcing at cubic order is 



N3 = J^Hr2o + ^ JlHr22 + ^ J^\HfH/2 + ^ Jl{AfHA + ArHA + Btha 
+BfHB + AfHB + BrHB))e'''^*-''^ 

+{^"jkJoAr2o + ^" JkJ2kAr2^ + ^" JkJoBr^^ + ^" JkJ2kBr^^ + ^"'\AfA/2 

+ ^"JoMHtha + Hthb + HrHA))e^-'^'' 
+($" Jfc Jo^r^0 + '^"JkJ2kAr^^ + <^"jkJoBr2o + <^" JkJ2kSr2^ + ^"'\B\'^B/2 

+^"'\A\^B + ^"joJkiHrHB + Htha + HrHB)e^+"^'') 
+ ... 
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Eliminating terms with dependencies e'"*, e"^ and e*"^ yields the three coupled amplitude equations, Eqs 28a||28c 



OtH = {n + in)AJoH + 2{a + iP)[C- 



\A\'+\BnH+HAB], 



OtA = {^i + m)/^JuA+{a + ib)\A\^A+{c+id)\B\^A+{a + il3)[2\H\^A + H'^B], 
OtB = {fi-in)AJkB+{a-ib)\B\^B + {c-id)\A\^B + {a-i^)[2\H\^B + H'^A], 

where a + if3, a + ib and c + id are given by Eqs|A5| |A6| and |A7| respectively. 



In the small delay limit (D 0) we can use the asymptotic values given by Eqs 12 to obtain, to leading order. 



a + ib 



id 




{D^)^[ 5 2 1-$'J2fc 



7r<i> 



+ i 



where x = '''■^/(^(l + ^) and a + i/3 is given by Eq, 
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(<&')^(6 + 7r) J2fc(<i> )- 

5 l-«''J2fc 



(A8) 
(A9) 



r 



Solutions and their stability 

Oscillatory uniform OU: 

This solution can be expressed as {H, A, B) 
(He*'^*,0,0), where 



n 



AJn 



The stability of this solution can be studied with the 
ansatz 

{H,A,B) = {'He''^\l + 5H+e^' + 5H^e^'), 
e-*'^*((5B+e^* + 5B_e^*)), 



which leads to three pairs of coupled linear equations 
which determine the six eigenvalues A. The first pair is 
restricted to the linear subspace of the small amplitude 
limit cycle and results in the standard stability problem 
which yields one stable eigenvalue A — — /iAJo and one 
zero eigenvalue corresponding to a shift in the phase of 
the oscillation. The other two pairs, which span the sub- 
spaces of ((5^+, (55+) and {5A-,5B-) respectively, give 



/ A-/xAJi-2aH2 + i(r2(AJo-AJi)-/37^2) -(^a + if3)n^ \fSA+\_ 

\ -{a - if3)H'^ A - /lAJi - 2an^ ~ i{n{AJo - AJi) - (3n^) ) \6B+ ) 

and the complex conjugate matrix spanning (SA^ , 6B^ ) . Setting the determinant equal to zero yields the characteristic 
equation 

A^ - 2A(A Ji - 2A Jo)m + A*^(A Ji - 4A JiA Jq + 3A Jq) + n^{AJi ~ AJ^f ~ 2-^VlAJq{AJ^ ~ A Jq) = 0. 

a. 

I 

There is an oscillatory instability for AJ\ = 2A J2 while This solution can be expressed as (iJ, A, B) = 

a steady instability occurs for AJ\ — AJq. The steady 
instability therefore always precedes the oscillatory one. 

Traveling waves TW: 
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(0,^Twe'"*,0), where 



A 



TW 



= n 



a / 

The stabihty of this solution can be studied with the 
ansatz 

{H,A,B) = (e"^*((5i?+e^* + (5iJ_e^*), 

which results in four coupled linear equations correspond- 
ing to the stability problem for TW in the competition 
between SW and TW (see section D, Turing-Hopf Bi- 
furcation). Here we assume that the TW solution is su- 
percritical and stable. We then turn our attention to the 
remaining two linear equations which describe the growth 
of the oscillatory uniform mode. These equations are un- 
coupled and yield the complex conjugate eigenvalues 

A = -A*(2^AJi-AJo)±*(f^(l-^) + (6~2/3)£)AJo, 

from which it is easy to see that an instability occurs 



for AJi = ^A Jq. This instability will generically occur 
with non-zero frequency. 

Standing waves SW: 

This solution can be expressed as {H, A, B) = 
(0, y^svye*'^', ^SH'e"*'^'), where Asw and w are given by 



equations (Eqs 31 and 32 1 




The stability of this solution can be studied with the 
ansatz 

{H,A,B) = {e"^\SH+e^* + 5H^e^'), 

Aswe'^\l + 6A+e^' + M_e^*), 



A 



This ansatz results in four coupled equations for the sta- 
bility of SW in the competition between SW and TW. 
Here we assume that the SW solution is supercritical and 
stable. The remaining two equations describe the growth 
of the oscillatory uniform mode. 



f X + iuj-{^i + in)AJo-AAsw{a + i(3) -2Asw{a + if3) \fSH+\ 

\ -2Asw{a-i(3) X-iuj-{fi-in)AJo-4:Asw{a-i(3) J \SH_ J 

Setting the determinant to zero yields the characteristic equation for the eigenvalues 

A2 - 2^A(A Jo - 4-^A Ji) + ^2(AJo - A-^AJ,f + (a Ji(r! - + ^^m) - ^AJof 

(a + c) {a + c) \ (a + c) [a + c) I 

^ i^AJl^a' ^(?)^^. (All) 



(a + c) 



The conditions for oscillatory and steady instabilities, Eqs [33] and |34) are found by setting A equal to iio and 
respectively. 

Mixed Mode: 

Mixed mode solutions are found by applying the ansatz Eq 36 to Eqs 28ap8c This gives 

n = nAJon + a{n^ + 2A^ + 2B^)n + 2nAB{acos(j)- Psin(j)), 
A = nAJiA+aA^ + cB^ A + 2an'^ A + n'^B{a cos (f> + 13 sin (f>), 
B = nAJiB + aB^ + cA'^B + 2an^B + n^A{a cos (I) + P sin (I)), 

= 2n{AJi-AJo) + 2f3n^ + {b+d-4p){A''+B^)-asin<j)(^^^ + ^^+AAB'^ 

+P cos cj^[^ + ^-4AB), (A12) 

where (j) = ipA — ipB ~ 26. One steady state solution of these equations takes the form {H, A, B, (j>) = (7i, A, — A, (/)), 
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where 



-[i^Jf){a + c) — /iA Ji ^ — 4q; + 2(q!C0S — /3sin0)^ 



a(a + c) — 


4a — 2(a cos 4> — fi sin 0) 
jA Jo(— 2q! + (a cos (p — (3i 


2a — {a cos (f)~ (3 sin </>) 
jin (/))) — /iA Jio; 


a(a + c) — 


4a — 2{a cos cj) — P sin 0) 


2q! — (a cos <j) ~ (3 sin </>) 



(A13) 
(A14) 



and (j) is found by plugging Eqs 



A13 



and 



A14 



study the stabihty of the mixed-mode state here. 



into Eq A12 and setting the left hand side equal to zero. We do not 



b. Double zero eigenvalue: Turing, Hopf 

z?(e): The solutions to the linear equation are time periodic oscillations and spatially periodic functions 

ri = H{T) + e''^* + A(T)e*'^'"= + c.c. 
i?(e^): The nonlinear forcing and resulting second order solution are 

N2 = J2(ff2e2^'^(*--f') + c.c. + 2\H\^)/2 + JliA^e^''''' + c.c. + 2\A\'^)/2 
<i>" JoJkHAc''=+"^' + ^" JaJkAHe-"^'+''= + c.c, 



rAH 



' AH 



r20 



-icj + l- ^'JkC^'^^ 



-"^^AH, 



'AH, 



<f"iJi\H\^+jim 



1 - Jo*' 



')9(e^): The nonlinear forcing at cubic order is 

= (^^"{J^Hr^Q + JlHr^H + + JIAtah) + ^"'4Jk\H\^A + <^"' JI\A\'' AI2^ e 

+ ( ^"{J^JkHrAH + JoJkHrAH + J^JkAr^o + JkJ2kAr2A) + J^J^H^^ A 



+<^"' Jl\A\^Al2\e''- 



Eliminating terms with dependencies e"^*, c^^ yields the two coupled amplitude equations, Eqs 39a|39b 

OtH = {fi + in)AJoH + {a + iP)\H\^H + {K + iA)\A\^H, 
OtA = fi^JkA + T\A\'^A + a\H\'^A, 

where /j, + ifi, a + ifi, fj and T are given by Eqs |A4| |A5[ |18a| and |18b] respectively, and 



-iujD 



K + iA 



1 + D^'Jqc-"^^ \ 1 - Jo Jie-^'^^ 



a> Jo J2 , 



1 



l + L'^'Ji 



— iujD 



-iw + 1 - Jie^"-° ' 2CJ + 1 - $' Jie-'^-oy ' l-$'Jo 



(A15) 



~*"'joJi 



(A16) 
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In the small delay limit {D — > 0) we can use the asymptotic values Eqs 12 to obtain, to leading order, 

TT {n/2 + i) ( - 



2£)($')3 (1 + 7r74) \ 



1 / J2fc($-) 



r = 



2 2(1-J2fe*')^ 

/ (ci>")2 ..A 



d{D), 



4i:>2($')3 \ 




and /i + ifi and a 
respectively. 



i/3 are given by Eqs 20a and 20b 



1e+06 



100 
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FIG. 15: Top: The real part of the cubic coefficient at the 
codimension-2 point. The sohd fine is the full expression, 
Eq |A6| and the dotted line is the asymptotic result in the 
Z) — > limit, Eq |A8| Bottom: The real part of the cross- 
coupling coefficient at the codimension-2 point. The solid 
line in the full expression, Eq |A7| and the dotted line is the 
asymptotic results in the D — > limit, Eq |A9[ 
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